function constraint_index = rome_conditional_constraint(cstr, varargin)
% ROME_VAR\ROME_CONDITIONAL_CONSTRAINT Adds a conditional constraint into
% the current model
%
%   constraint_index = rome_conditional_constraint(cstr, Map1, Map2, ...)
%
% Input:
%     cstr: Constraint. Must be an LDR inequality (in the nonnegative
%     orthant cone).
%     Map1, Map2, ... : Map(s) used to clone uncertain variables.

% check the pathological case
if (isempty(varargin))
    constraint_index = rome_constraint(cstr);
    return
end

% check input validity
if (~cstr.IsLDR || ~isscalar(cstr) || cstr.Cone ~= rome_constants.NNOC)
    error('rome_var:rome_conditional_constraint:InvalidConstraint', ...
        'Constraint must be a scalar LDR inequality!');
end

% obtain the primitive index set for cstr
cstrMap = reshape(cstr.BiAffineMap, cstr.NumMappedVars + 1, ...
    cstr.NumMappedRandVars + 1);
cstrPrimInd = cstr.NumUnmappedRandVars + ...
    find(any(cstrMap(:, 2:end), 1));

% create the primitive index set for the conditional constraints
zPrimInd = cstrPrimInd;
for ii = 1:length(varargin)
    Map = varargin{ii};
    zPrimInd = [zPrimInd, Map.origPrimInd];
end
zPrimInd = unique(zPrimInd);

% get the current model
h = rome_get_current_model();

% obtain rand vars
z = [1; h.get_rand_vars(zPrimInd)];

% create dual variables for cstr
r = rome_model_var(length(zPrimInd)+1);

% construct dual constraints for cstr
dual_sum = r;
dual_cstr = r'*z;

% construct dual constraints for conditional cstrs
for ii = 1:length(varargin)
    Map = varargin{ii};
    
    % make additional rand var clones
    AddPrimInd = setdiff(zPrimInd, Map.origPrimInd);
    if ~isempty(AddPrimInd)
        [zAdd MapAdd] = rome_clone_uncertainty(AddPrimInd);
        
        [tmp id] = sort([Map.origPrimInd, MapAdd.origPrimInd]);
        zzPrimInd = [Map.clonedPrimInd, MapAdd.clonedPrimInd];
        zzPrimInd = zzPrimInd(id);    
    else
        zzPrimInd = Map.clonedPrimInd;        
    end
    
    zz = [1; h.get_rand_vars(zzPrimInd)];
    % create dual variables for cstr
    s = rome_model_var(length(zPrimInd)+1);
    
    % construct dual constraints for cstr
    dual_sum = dual_sum + s;
    dual_cstr = dual_cstr + s'*zz; 
end

% obtain the LDR coefficients for cstr
x = h.get_vars(cstr.NumUnmappedVars + (1:cstr.NumMappedVars));
ycstr = cstrMap' * [1; x];

% transform the coefficients for cstr to that for z
y = rome_var(length(zPrimInd)+1);
y.NumUnmappedVars = ycstr.NumUnmappedVars;
y.BiAffineMap = sparse(length(zPrimInd)+1, ycstr.NumMappedVars+1);
IndInCstrFlag = (zPrimInd > cstr.NumUnmappedRandVars) & ...
    (zPrimInd <= cstr.NumUnmappedRandVars + cstr.NumMappedRandVars);
y.BiAffineMap([true IndInCstrFlag], :) = ycstr.BiAffineMap( ...
    [1, zPrimInd(IndInCstrFlag) - cstr.NumUnmappedRandVars + 1], :);

% impose dual constraints
rome_constraint(dual_sum == y);
constraint_index = rome_constraint(dual_cstr >= 0);